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ABSTRACT 

Instantaneous  frequency  (IF)  estimation  of  signals  with  non¬ 
linear  phase  is  challenging,  especially  for  online  processing. 
In  this  paper,  we  propose  IF  estimation  using  sequential  Baye¬ 
sian  techniques,  by  combining  the  particle  filtering  method 
with  the  Markov  chain  Monte  Carlo  (MCMC)  method.  Us¬ 
ing  this  approach,  a  nonlinear  IF  of  unknown  closed  form 
is  approximated  as  a  linear  combination  of  the  IFs  of  non¬ 
overlapping  waveforms  with  polynomial  phase.  Simultane¬ 
ously  applying  parameter  estimation  and  model  selection,  the 
new  technique  is  extended  to  the  IF  estimation  of  multicom¬ 
ponent  signals.  Using  simulations,  the  performance  of  this 
sequential  MCMC  approach  is  demonstrated  and  compared 
with  an  existing  IF  estimation  technique  using  the  Wigner  dis¬ 
tribution. 

Index  Terms —  Frequency  estimation,  Bayes  theorem,  par¬ 
ticle  filter,  Markov  chain  Monte  Carlo 

1.  INTRODUCTION 

Many  naturally  occurring  signals  have  time-varying  (TV)  spec¬ 
tral  characteristics.  Whereas  sinusoids  have  a  single  frequency 
over  all  time,  linear  frequency-modulated  (FM)  chirps  (used 
in  radar  and  sonar)  have  frequencies  that  vary  linearly  with 
time,  and  shallow-water  acoustic  waves  undergo  dispersive 
(nonlinear)  frequency  modulations. 

For  a  single  component  signal, 

s(t)  =  A{t)e^cv>w,  (1) 

with  amplitude  A(t),  phase  p{t )  and  FM  rate  c,  the  frequency 
at  a  particular  time  can  be  described  by  the  IF: 

C s{t)  =  C^(f).  (2) 

In  many  real-life  applications  such  as  radar,  sonar,  under¬ 
water  acoustics  and  structural  health  monitoring,  the  IF  can 
be  a  powerful  tool  in  estimating  important  parameters  of  the 
signal.  For  example,  it  was  shown  in  [1]  that  varying  C,s(t) 
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from  a  linear  to  a  hyperbolic  function  resulted  in  better  track¬ 
ing  estimation.  Furthermore,  the  ability  to  estimate  the  IF  of 
a  signal  propagating  through  an  unknown  environment  may 
provide  information  useful  in  characterizing  the  environment. 
As  a  result,  it  is  important  to  be  able  to  accurately  estimate 
the  IF  to  obtain  the  signal’s  spectral  variation  with  time. 

In  this  paper,  we  begin  with  a  brief  overview  of  the  cur¬ 
rent  IF  estimation  methods  and  their  shortcomings  in  Section 
2.  In  Section  3,  we  propose  a  new  IF  estimation  method  based 
on  the  use  of  sequential  Bayesian  techniques.  Specifically,  we 
combine  a  particle  filter  with  MCMC  to  estimate  static  param¬ 
eters  of  the  IF  of  windowed  signal  segments.  We  extend  this 
to  signals  with  multiple  components  using  model  selection. 
Some  simulation  results  that  demonstrate  the  effectiveness  of 
the  new  approach  are  presented  in  Section  4. 

2.  IF  ESTIMATION  METHODS 

Parametric  and  non-parametric  approaches  to  IF  estimation 
are  commonplace  in  the  literature.  Non-parametric  IF  esti¬ 
mation  techniques  do  not  assume  a  mathematical  model  for  a 
signal  and  instead  use  time-frequency  (TF)  signal  representa¬ 
tions  to  describe  the  IF  in  the  TF  plane.  Conversely,  paramet¬ 
ric  IF  estimation  methods  assume  a  specific  signal  model,  as 
in  (1),  and  then  proceed  to  estimate  the  IF. 

TF  analysis  is  a  well-known  processing  tool  that  has  been 
applied  to  the  IF  estimation  of  signals  as  a  non-parametric 
approach  as  it  provides  a  natural  means  to  display  the  sig¬ 
nal  spectrum  at  every  time  instant.  The  Wigner  distribution 
(WD),  Ws(t,  f )  =  /  s(t  +  r/2)s*(t  —  t /2)e-i27rfT dr,  of  a 
signal  s{t)  is  a  very  popular  time-frequency  representation 
(TFR)  that  can  provide  high  resolution  along  the  signal’s  IF. 
In  particular,  it  can  be  shown  that  the  IF  in  (2)  of  the  signal  in 
(1)  is: 

C»(t)  =  J  f  Ws(t,f)df.  (3) 

For  multicomponent  signals,  Ai{t)e^2vc(piY\  with  nonlin¬ 
ear  phase  i pi{t),  the  IF  in  (2)  no  longer  provides  the  signal’s 
matched  TV  spectrum  since  (2)  is  not  a  linear  representation. 
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This  follows  also  from  (3)  since  the  WD  is  a  quadratic  TFR 
that  suffers  from  cross  terms  for  multicomponent  signals. 

The  IF  can  also  be  estimated  by  extracting  ridges  (or  peaks) 
of  a  TFR  and  then  applying  a  peak  detection  algorithm  to  es¬ 
timate  the  IF: 

C s(t)  =  argjmax  {Ts(t,  /;  /i)}}, 

J 

where  Ts  (t,  /;  h)  is  a  quadratic  TFR  that  depends  on  a  smooth¬ 
ing  window  h(t)  and  T  contains  frequency  values  for  which 
the  TFR  is  non-zero.  An  algorithm  that  used  a  smoothed  WD 
with  an  adaptive  window  width  determination  to  improve  the 
estimation  performance  was  used  in  [2],  Although  this  can 
improve  the  performance  of  the  WD  for  multicomponent  sig¬ 
nals,  the  smoothing  may  smear  the  estimated  IF  values.  IF 
estimation  of  seismic  data  was  used  in  conjunction  with  the 
reassignment  method  in  [3],  A  drawback  of  the  original  reass- 
ingment  method  was  a  loss  of  performance  in  low  signal-to- 
noise  ratio  (SNR)  environments.  To  mitigate  noise-sensitivity, 
a  multi -taper  TF  reassignment  method  was  proposed  in  [4]. 

IF  estimation  methods  have  also  been  developed  based  on 
parametric  and  often  statistical  models  to  describe  the  sig¬ 
nal.  For  example,  Newton’s  method  was  used  in  [5]  to  de¬ 
termine  the  maximum  of  the  log-likelihood  function,  but  this 
only  worked  for  the  joint  estimation  of  frequency  and  FM  rate 
of  a  signal  in  white  Gaussian  noise.  In  [6],  a  Bayesian  ap¬ 
proach  for  IF  estimation  was  proposed  by  using  the  MCMC 
method  with  parametric  processes  when  prior  knowledge  on 
the  statistical  properties  of  the  processes  was  available. 

3.  SEQUENTIAL  MCMC  IF  ESTIMATION 

The  previously  mentioned  methods  do  not  perform  well  when 
(i)  the  IF  is  nonlinear,  (ii)  a  closed  form  expression  for  the  sig¬ 
nal  is  unavailable,  or  (iii)  the  signal  is  multicomponent.  More¬ 
over,  the  estimates  may  need  to  be  computed  offline;  in  a  wide 
variety  of  applications,  estimates  are  required  to  be  updated 
online  and  executed  sequentially.  Finally,  the  implementa¬ 
tion  method  may  require  storage  of  intermediate  data  leading 
to  excessive  storage  costs  and  increased  computational  com¬ 
plexity  as  additional  measurements  are  computed. 

In  this  work,  we  propose  a  sequential  Bayesian  IF  estima¬ 
tor  that  combines  a  particle  filter  with  MCMC  methods  for  a 
non-parametric  signal  model.  Without  assuming  the  existence 
of  a  closed  form  expression  for  the  phase  <p(t),  we  window 
the  signal  into  non-overlapping  segments  with  each  segment 
modeled  with  polynomial  phase  with  unknown  coefficients. 
Thus,  we  approximate  the  IF  ( s(t)  as  a  linear  combination  of 
the  derivatives  of  the  polynomial  functions.  The  static  coeffi¬ 
cients  are  estimated  sequentially  for  each  window:  a  particle 
filter  allows  for  online  processing  and  the  MCMC  increases 
estimation  accuracy  and  reduces  storage  cost  and  computa¬ 
tional  complexity.  Since  the  IF  in  each  window  can  be  ap¬ 
proximated  by  a  polynomial,  the  problem  is  reduced  to  deter¬ 
mining  the  polynomial  coefficients.  As  expected,  the  larger 


the  degree,  the  better  the  estimation  accuracy.  However,  there 
is  a  trade-off  between  accuracy  and  computational  complex¬ 
ity. 

3.1.  Estimation  of  Static  Parameters 

As  particle  filtering  was  designed  for  dynamic  parameter  es¬ 
timation,  it  was  shown  to  fail  when  used  to  estimate  static 
parameters  [7].  Following  the  method  discussed  in  [8],  we 
estimate  the  parameter  vector  in  each  window  by  combining 
sequential  importance  sampling  particle  filtering  with  MCMC 
methods.  Next,  we  provide  the  main  steps  of  this  sequential 
MCMC  (SMCMC)  approach. 

Suppose  at  time  step  k,  the  particles  and  their  correspond¬ 
ing  weights,  (x1,  wk),  (x2,  to2),  •  •  •  ,  (xNe,  w^3),  are  used  to 
represent  the  conditional  probability  density  function  p{x\Zk), 
where  Ns  is  the  number  of  particles,  x  is  the  static  parameter 
vector  to  be  estimated,  z k  is  the  observation  vector  at  time 
step  k  and  Zk  =  [zi,  Z2,  •  •  •  .  zk],  The  SMCMC  method  up¬ 
dates  the  weight  for  each  particle  using: 

wk+m  ocp(zfc+i,zfe+2,--  -  ,zk+m\xl,Zk)wl,  i  =  1,  •  •  •  ,NS 

where  m  is  the  batch  size.  Next,  a  rejuvenation  test  is  per¬ 
formed  using  the  Kullback-Leibler  distance  measure  [8].  Spe¬ 
cifically,  if  severe  degeneracy  occurs  [8],  indicated  by  the  test, 
the  MCMC  method  is  applied  using  the  independent  Metropo- 
lis-Hastings  (IMH)  algorithm  [8]  with  the  Gaussian  distribu¬ 
tion  A/"(x|/ox,  £x)  as  the  proposal  density,  where 

Ne  Ns 

Mx  =  ^  '  Wk+mX  i  =  ^  1  Wk+m(X  —  Afx)(x  —  Mx)  • 

i—1  !-] 

New  particles  and  weights  are  then  obtained  by  sampling  from 
this  Gaussian  distribution,  representing  p(x.\Zk+m). 

When  x  can  originate  from  different  types  of  signal  struc¬ 
tures  {Hi,H2,  •  •  •  ,  Hm},  model  selection  is  used  simultane¬ 
ously  [8].  Let  the  parameter  vector  for  model  j  be  then, 
the  resulting  distribution  can  be  obtained  using: 

M 

p(x|Zfc)  =^P(Hj\Zk)p(xW\Zk,Hj),  (4) 

j= i 

where  P(Hj\Zk)  is  the  probability  of  model  Hj  given  mea¬ 
surements  Zfc,  and  p(xfj) \Zk,  Hj)  is  the  one  discussed  above 
given  model  Hj .  Instead  of  using  multiple  hypothesis  testing 
in  [9],  the  model  probabilities  are  updated  sequentially. 

3.2.  IF  Estimation  of  Single  Component  Signals 

The  IF  of  a  single  component  signal  (that  can  be  decomposed 
into  basis  functions  which  are  non-overlapping  in  time)  is 
given  by  the  linear  combination  of  the  IFs  of  these  basis  func¬ 
tions.  Specifically,  assume  that  s(t)  can  be  decomposed  as: 

L 

s(t) 

1  =  1 


where  pi{t)  =  u(t  —  (/  —  1)T)  —  u(t  —  IT)  is  a  rectangular 
window  and  u(t)  is  the  unit  step.  Then  we  can  show  that  the 
IF  of  s(t)  is  given  by  C s(t)  =  E^=  1  There¬ 

fore,  given  a  waveform  with  an  unknown  nonlinear  phase 
we  propose  to  approximate  its  IF  £s(t)  using  a  linear 
combination  of  IFs  of  FM  waveforms  with  polynomial  phase 
tpi(t)  =  Elo  an,itn  and  duration  T.  Specifically, 

L  N 

Cs(t)  ~  '^2'^2nantitn~1pi(t).  (5) 

1=1  n=  1 

With  this  assumption,  the  estimation  of  the  IF  becomes  the 
estimation  of  a  set  of  unknown  static  parameter  vectors  x/  = 
[at, z,  0,2,1,  ■  ■  •  ,  ajv.z],  l  =  1,  •  ■  • ,  L  using  the  SMCMC  method 
described  in  Section  3.1.  In  each  window  /,  a  new  set  of  parti¬ 
cles  is  used  to  estimate  the  parameter  vector  x(.  To  reduce 
approximation  errors,  we  consider  short  duration  windows 
pi(t)’,  the  window  length,  however,  has  to  be  long  enough  to 
reduce  errors  in  estimating  the  polynomial  coefficients. 

3.3.  IF  Estimation  of  Multiple  Component  Signals 

For  a  multicomponent  signal,  the  IF  cannot  be  simply  inter¬ 
preted  as  the  average  frequency  at  each  time  [10].  Instead,  we 
need  to  obtain  a  representation  that  provides  the  sum  of  the 
IFs  of  the  different  signal  components  in  order  to  provide  the 
frequency  content  of  the  signal.  If  we  assume  that  the  num¬ 
ber,  Q,  of  components  of  a  multicomponent  signal  is  known, 
and  that  the  signal  can  be  decomposed  as: 

q=  1  1=1 

where  A (t)  and  (t )  are  the  amplitude  and  phase  of  the 

gth  component  in  the  /th  window,  respectively,  then  the  IF  of 
this  component  is  approximately: 

N 

(6) 

n=  1 

where  a^\  is  the  ?rth  coefficient  of  component  q  in  the  /th 
window.  Therefore,  in  the  /th  window,  Q  parameter  vectors, 
xj9^  =  [a1]9; ,  ■  ■  •  ,  a^],  q  =  1,  •  •  ■  ,  Q,  are  estimated  using 
SMCMC  with  a  new  set  of  particles  for  each  vector. 

Our  approach  can  be  extended  to  the  IF  estimation  of  mul¬ 
tiple  component  signals  when  the  number  of  components  is 
unknown.  In  this  case,  we  assume  that  there  are  several  mod¬ 
els  corresponding  to  the  number  of  components  present  in 
the  /th  window.  For  each  model,  Qi  FM  waveforms  (Qi  £ 
[1,  2,  •  •  •  ,  Q])  are  used  to  approximate  the  IFs  of  the  Qi  com¬ 
ponents  in  the  /th  window  with  Qi  sets  of  polynomial  coef¬ 
ficients  x|9^  =  [a[9\  •  ’  ’  9  =  1,  •  •  ■  ,  Qi  as  the  pa¬ 

rameter  vectors  to  be  estimated.  The  number  of  components 


Qi,  which  may  vary  in  different  windows,  is  determined  us¬ 
ing  model  selection.  Specifically,  if  the  maximum  number  of 
components  is  known  to  be  Q ,  there  are  totally  Q  +  1  models 


to  choose  from: 

H0: 

zi(t) 

=  Vi(t) 

Hi  : 

Zl  (t) 

=  +  Vl{t) 

Hq: 

=  £?=i  A^ity^w  +  Vl(t) 

The  Qi  signal  component  parameter  vectors  in  the  /th  window 
are  estimated  simultaneously. 

4.  SIMULATIONS 

With  the  order  N  =  2  in  (5),  linear  approximation  works  best 
when  there  are  no  sudden  changes  in  the  IF  and  if  the  window 
length  is  small.  We  demonstrate  the  estimation  performance 
by  comparing  the  true  and  estimated  IF  of  a  9  dB  SNR  noisy 
single  component  signal  in  Fig.  1  (a)  with  window  length  500 
samples.  The  time  averaged  root  mean  square  error  (RMSE) 
for  varying  SNR  is  shown  in  Fig.  1  (b). 


(a)  (b) 

Fig.  1.  IF  estimation  using  linear  approximation:  (a)  esti¬ 
mated  compared  with  true  IF,  (b)  RMSE  vs  SNR. 

Quadratic  approximation  (polynomial  with  order  N  =  3) 
is  performed  in  Fig.  2  (a),  and  compared  with  linear  approx¬ 
imation.  We  used  a  9  dB  SNR  signal  and  1000  samples  win¬ 
dow  length.  It  can  be  seen  that  for  this  case,  where  there  is  a 
sudden  change  in  the  IF,  the  quadratic  approximation  is  much 
better  than  the  linear  one,  resulting  in  a  smaller  estimation 
error.  The  averaged  RMSE  is  4.13  Hz  for  the  linear  approxi¬ 
mation  and  0.80  Hz  for  the  quadratic.  Note  that  the  window 
length  was  chosen  large  in  this  example  for  comparison.  The 
larger  the  window  length,  the  worse  the  linear  approximation. 

Next,  we  demonstrate  the  use  of  SMCMC  with  model  se¬ 
lection  for  multicomponent  signal  IF  estimation.  Fig.  3  shows 
the  estimation  of  two  component  signals  overlapping  in  time 
as  they  approach  each  other  in  frequency.  The  SNR  is  9  dB 
and  the  window  length  is  500.  The  performance  demonstrates 
the  advantage  of  our  method  even  when  the  two  components 


5.  CONCLUSION 


Fig.  2.  (a)  Comparison  of  IF  estimation  of  a  single  com¬ 
ponent  signal  using  quadratic  and  linear  approximation,  (b) 
IF  estimation  of  two  crossing  linear  FM  chirps  using  linear 
approximation  SMCMC  (stars  and  crosses),  compared  with 
short-time  WD  (dots). 


are  very  close.  Note  that  the  IF  estimation  in  Fig.  3  (b)  is  very 
difficult  for  most  IF  estimation  methods. 


(a) 


Fig.  3.  SMCMC  IF  estimation;  signal  components  are  (a)  not 
overlapping  and  (b)  overlapping  in  frequency. 


We  compare  the  performance  of  our  approach  with  the 
WD  IF  estimation  technique.  For  comparison,  the  method 
in  [11]  is  extended  to  short-time  processing  by  executing  the 
WD  IF  estimation  on  windowed  data  segments.  Specifically, 
we  estimate  linear  FM  parameters  of  the  WD  of  windowed 
segments  that  are  512  samples  long  using  the  method  in  [11]. 
The  simulation  result  is  shown  in  Fig.  2  (b)  for  9  dB  SNR. 
The  WD  suffers  from  cross  terms  when  it  is  applied  to  mul¬ 
tiple  component  signals.  If  two  curves  are  close  in  the  TF 
plane,  the  WD  IF  estimation  method  can  pick  the  cross  term 
location  as  the  signal  component;  this  is  not  the  case  with 
the  SMCMC.  Note  that  to  implement  the  WD  estimation,  the 
number  of  signal  components  has  to  be  known;  only  the  max¬ 
imum  number  of  components  is  assumed  known  for  the  SM¬ 
CMC.  Also,  the  computational  complexity  and  storage  cost 
are  reduced  dramatically  by  our  approach.  The  new  method 
takes  less  than  half  of  the  CPU  time  and  a  quarter  of  the  stor¬ 
age  cost  when  compared  to  the  short-time  WD  approach. 


In  this  paper,  we  proposed  an  approach  for  the  IF  estimation 
of  both  single  and  multiple  component  signals  with  nonlin¬ 
ear  phase  that  may  not  exist  in  closed  form.  Our  approach 
approximates  the  IF  of  a  signal  by  a  linear  combination  of 
the  IFs  of  windowed  FM  waveforms  with  polynomial  phase. 
Although  we  have  shown  results  comparing  our  method  to 
the  short-time  WD,  we  plan  to  consider  adaptive  windows  in 
the  SMCMC  (to  improve  performance  in  lower  SNRs)  and  to 
compare  our  method  with  the  multi-taper  reassignment  in  [4]. 
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